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Abstract 

The flight trajectory of aerospace vehicles subject to a class of path constraints is 
considered. The constrained dynamics is shown to be a natural two-time-scale system. 
Asymptotic analytical solutions are obtained. Problems of trajectory optimization and 
guidance can be dramatically simplified with these solutions. Applications in trajectory 
design for an aerospace plane strongly support the theoretical development. 


1. Introduction 

In many flight control and trajectory optimization problems, certain portions of the"' 
trajectory are required to follow some state space constraints dictated by operational or 
safety considerations. Nowadays almost all realistic atmospheric flight problems are numer- 
ically intensive. Solutions to these constrained problems are thus in the form of numerical 
data. Analytic expressions of the constrained trajectory as explicit functions of time, if 
available, are the dream of a trajectory designer. They provide an efficient means to eval- 
uate the trajectory, and often lead to a better understanding of the trajectory. In turn, 
tasks such as trajectory optimization, control and guidance can be significantly simplified. 
In some cases approximate analytical solutions may be possible because of the additional 
relationships due to the constraints. Gilbert et al present an enlightening treatment of 
a coasting arc observed in the optimal trajectory of a ground-based interceptor 1 . But in 
general, no systematic approach exists to obtain analytical solutions to the constrained 
system which is often nonlinear in nature. This paper extends the approach in Ref. 1 
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to more general flight problems. For an important class of constrained flight problems 
which axe described in Section 2.1, the dynamics is shown to be a two-time-scale system. 
Approximate analytical asymptotic solutions are obtained. Successful applications in the 
optimal ascent and hypersonic cruising trajectory analyses for an aerospace plane are pre- 
sented. Comparison of the approximate solutions with numerically generated solutions 
shows excellent agreement. 

2. Theoretical Development 

2.1 Problem Formulation 


Consider the point-mass model for a vehicle flying in a fixed great circle plane over a 
nonrotating spherical earth. Define the dimensionless variables 


h = 



V t-t 0 

\f9oT~o \J r o! 9o 


( 1 ) 


*^where r is the distance from the center of the earth to the vehicle; V the speed; t the 
^current time. r 0 and g 0 are the radius and gravitational acceleration, respectively, at the 
starting time t 0 . The dynamics is given by 


h' = v sin 7 
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The prime in Eqs. (2)-(5) denotes the differentiation with respect to r. An algebraic 
constraint is imposed on the trajectory 


C(v,p(h)) = 0 


( 6 ) 


In the above equations 7 stands for the flight path angle, 8 the down range polar angle, T 
the thrust, L and D the aerodynamic lift and drag, respectively, p in Eq. (6) represents 
the atmospheric density. 
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Assumptions : 

(i) The dependence of the magnitude of the thrust T on altitude, if any, is exclusively 
through atmospheric density or dynamic pressure. 

(ii) The atmospheric density is an exponential function of altitude 

p = p o e~ 0h = p a p(h ) (7) 


where j3 = r 0 /h a with h 3 being the scale height. 

(iii) The constraint contains p explicitly, i.e., dC/dp 7 ^ 0 . At the starting point 
C(v(t 0 ), p 0 ) = 0. dC/dv 7 ^ 0 for t > t 0 . Then by the implicit function theorem v can be 
regarded as a function of p 

v=nm) ( 8 ) 

(iv) The set of feasible controls is not empty. We define here a pair of controls 
[T(<), «■(<)] to be feasible if the solutions of Eqs. (2)-(5) satisfy Eq. ( 6 ), and 0 < T(t) < 
T m ax and a m j„ < a(t ) < o max , where T max is the maximum available thrust (which may 
be altitude dependent), and a m ,„ and a max are the lower and upper bound of angle of 
attack, respectively. 

(v) I 7 I is small along all feasible trajectories. 

Assumption (i) is valid for many flight scenarios and both airbreathing and rocket 
propulsion systems. Note that (i) includes pure coasting case when T = 0. Assumptions 
(ii) is easily justified. Assumption (iii) specifies the class of constraints under consideration. 
Many important physical constraints in hypersonic flight are in this class. Assumption (v) 
is often satisfied during constrained flight. 


2.2 Asymptotic Solutions 


In the following analysis, we assume that the vehicle is flying along a feasible trajectory. 
We proceed by first differentiating ( 8 ) once 


, df dp , . 

” =Tpdh h 


By (4) and (9), 


-Pfppv sin 7 = 


T — D sin~f 


mg 0 (1 + h) 2 


( 9 ) 


( 10 ) 
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There axe two controls for the vehicle, T and a. Either one of the controls can be pre- 
programmed or determined by some type of performance optimization. The other control 
is then determined from Eq. (10) to satisfy the constraint, provided that the resulting 
control pair is in the set of feasible controls defined by Assumption (iv). It is understood 
that the following discussions apply to such a case. Once T and o are determined, we can 
define the lift-to-drag ratio 


V = 


D 


and the thrust-to-drag ratio 


T 

k = d 


( 11 ) 


( 12 ) 


t] is primarily a function of Mach number for given a, and k a function of Mach number 
and atmospheric density for given throttle setting. Next, neglect the gravity component 
in the v' equation for small 7 


v = 


T-D (k-1 )D 


mg D 


mg 0 


(13) 


Consider the flight path angle equation (5). We note that typically h << 1 for atmospheric 
flight, and 7 is assumed to be small. Using cos 7 sa 1 , 1 + h ~ 1 and L = rjD in (5), we 
have 

1 n D 

(14) 


7 ' = I(_^__1 +u 2) 


v mg 0 

Assume that k ^ 1 . Replacing D in (14) by (13) yields 


/ 1 , v'n 2 


Noticing v = f(p), we substitute v 1 in (15) by Eq. (9) to obtain 
To explore the solution, we define new scaled variables 


(15) 


(16) 


h = fih, 7 = £7 


(17) 
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From Eqs. (2) and (16) the equations for the two new variables are 


( 18 ) 

(19) 


h' = vj 

7 *« 


where e = 1//3 = h 9 /r 0 . For the earth, e is on the order of 10~ 3 . By Assumption (i) k is 
a function of h through the density. The right hand side of Eq. (19) is independent of £. 
Therefore, we see that equations (18)— (19) constitute a natural two- time- scale system. By 
the standard singular perturbation theory 2 , the sufficient condition that fpT)/(K — 1) > 0 
guarantees that the system has an attracting asymptotic solution to which the solution of 
the system will quickly converge. The zeroth-order outer solution is obtained by setting 
the right hand of Eq. (19) equal to zero 


=, _*(/-/ 

7 = c t 

Jp *7 


( 20 ) 


where p = e h has been used. Substituting (20) and v = f into (18) gives 

dh 1 -/ 2 

— = -P J 
dr pfp 


where 


Noticing that df = — pfpdh , we have 



df 

1 -P 


= PdT 


which, upon integration, yields 


In (24) Q and z(r) are defined by 


e 2 ^ r > + Q 

e 2z(r) _ Q 

V 0 ~ 1 
V 0 + 1 

z(t ) = f P(a)da 
Jb 


V = f = 


Q = 


( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 
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Discussions 

(1) Equation (20) provides an approximate closed-form solution for 7 along the con- 
strained trajectory. This formula can be very useful for guidance purposes. The accuracy 
is to the degree of sin 7 w 7 and 1 + h « 1 . Notice that / and fp can be evaluated from C 
[Eq. ( 6 )] via the implicit function theorem. 

( 2 ) The solution (20) is the zeroth-order outer solution. In general, the initial condition 
on 7 will not be satisfied. In dimensionless time r, this solution will only be accurate 
outside the initial boundary layer of thickness on the order of e, which corresponds to a 
few seconds in real time. 

(3) When T = D (or k = 1) at some isolated points, a singularity arises at those points. 
Formulas (20) and (24) then represent the outer solutions in the intervals separated by 
those points. There is a small neigborhood around each of those points where the the outer 
solutions are inaccurate. When T = D, it can be shown from Eqs. ( 2 ) and (4) that 



If we substitute above equation into Eq ( 6 ), Eq. ( 6 ) becomes an algew*kic equation in v. 
It is clear that the constrained trajectory can only have a constant speed which is a root of 
Eq. ( 6 ). This in turn implies that the altitude can only be constant. Thus 7 = 0. Letting 
k = 1 and z(t) = 0 in Eqs. (20), ( 21 ) and (24), we arrive at 7 = 0, h = ho = 0 and v = v 0 , 
which verifies that the above formulas are still valid in this case. 

( 4 ) Interestingly, the variation in v is not explicitly dependent on the form of the 
constraint ( 6 ) according to Eq. (24). A similar formula is obtained in Ref. 1 for a specific 
constraint. But here we show that this is a general result for the class of constraints 
under consideration. Another nice feature of the above formulas is that there is no explicit 
dependence on the mass m which is time- varying for T ^ 0. The only inconvenience in 
the explicit formulas (24) is the function z(t). In some case z may take a simple form as 
an example will demonstrate later. For more general situations often an estimate of the 
average value of P, say, P, may be used in (26). Then, 

z(t) « Pr (27) 
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If necessary, P can always be estimated without much difficulty from numerical simulation 
of the accurate trajectory. The errors introduced by this approximation are limited because 
any errors in h and 7 are about one thousand times smaller than those in h and 7 , owing 
to the scaling effect of Eq. (17). In the post-run data analysis, having to estimate P in 
the analytic representations is much better than to sort through the numerical data. With 
this simplification the explicit solution for 9 is found to be 

6' & v =*► 0 - 8 0 = - r (28) 

( 5 ) Up to this point, the discussion is general in the sense that no specific forms of 
C in Eq. ( 6 ) have been assumed, as long as Assumption (iii) is valid, h is an implicit 
function of time through ( 6 ) and (24). If, however, h can be solved explicitly in terms of 
v from Eq. ( 6 ), then h can be expressed as a function of time. 7 can thus be obtained as 
an explicit function of time through ( 20 ). 

Let us consider two types of constraints that will be used in the next section. 

Case 1 

C = vp p — C p = 0, Cp > 0 and p any nonzero real number (29) 


Hence 

v = f = Cpp~ p = C p e p ~\ h=- In {f l Cp) 

P 

The solutions for h and 7 are 



e 2 ^ + Q 
C p (e^)-Q) } 



e 2 *< r > - Q 

e 2z(r) + Q 


) 2 -l) 


Case 2 

C is a function of v and p, and v can be expressed as 


(30) 


(31) 

(32) 


v = f = (A + Bp m ) n 


(33) 
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where A, B, m and n Eire all nonzero real numbers such that f(e~ h ) > 0. Case 1 can 
actually be included in (33) if we allow A to be zero. Nevertheless, we list the formulas 
here seperately for the convenience of later reference. 


r 


(e 2z < r > + Q) l / n - A(e 2z < r > - Q)^ n 
B(e 2dr) _ Q) l/n 


= 9(r) 


(34) 


7 (r) = P 


(A + Bg(r))” - (A ± Bg(r))~" 
nmB(A + B g(T)) n ~ l g(r) 


1 


h(r) = lng(r) 


m 


(35) 

(36) 


We conclude this section by stressing that the point at to ( T = 0) does not have to be 
the point where the trajectory just enters the constraint (6). It can be a midway point on 
the constraint. By letting r have negative values, the expressions for t/(r), 7 (r), h(r) and 
9(t) are valid for the portion of the trajectory before to- This feature can be very useful 
in characterizing a family of constrained trajectories which enter the same constraint at 
different points. 


3. Applications in Trajectory Analysis for an Aerospace Plane 
3.1 Optimal Ascent to Orbit 

An aerospace plane is a hypersonic vehicle that has single-stage-to-orbit (SSTO) ca- 
pability. The primary propulsion comes from airbreathing engines. There has been active 
study of the fuel-optimal ascent trajectory for the aerospace plane 3-6 . It is found that a 
typical optimal trajectory will sequentially climb on two operational constraints 


g — g 2 — Qmax 

(36) 

Q = Ky/pV 3 = Qmax 

(37) 


The first one is a dynamic pressure constraint and the second a convective heating rate 
constraint. In conformity with Eq. (29), the two constraints can be rewritten as, respec- 
tively, 


v? 5 = ve-° sh = C, 

(38) 

vp 1/6 =ve- h/< ‘ = C Q 

(39) 
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where 


Cq — \/^QmaxJ~po9o^o -i Cq — {Q max / P ffo^o^Poffo^o) 


(40) 


A generic aerospace plane model 7 is used. The thrust of the airbreathing propulsion 
system is modeled by 

T = C r q (41) 

where the thrust coefficient Cj is a function of dynamic pressure, Mach number and the 
throttle, controlled by the fuel equivalence ratio. Hence the model satisfies Assumption 
(i). For the best fuel efficiency the equivalence ratio is set at unity. Assumption (v) has 
been verified 6 . 

A typical optimal ascent trajectory for the aerospace plane is shown in Fig. 1. After a 
quick initial climb-out, the trajectory enters the ^-constraint boundary and rides on it until 
the Q-constraint becomes active. After flying on the Q-constraint awhile the aerospace 
plane finally pulls up at an appropriate point and ascends to orbit. When the aerospace 
plane is flying on the constraint (36), the complete trajectory is characterized by 


e 2 ^ + Q 

V ~ ~ 

e 2p ' - Q 

(42) 

, e 2pr + <? , 

(43) 

p 2Pt — O 

^2Pn elPr+ Q Q y-D 

(44) 

2Pr _ Q 

0 = e„ + p in( y _ q'-)~ r 

(45) 


The approximation of (27) has been used in the above equations. Similar formulas with 
p = 1/6 in Eqs. (31) and (32) describe the motion analytically after the aerospace plane 
leaves the constraint (36) and enters the constraint (37). When q ma x — 95,660 N/m 2 
(2,000 psf ) and Qmax = 800 W/cm 2 , we choose P = 0.9 and P = 1.2 for the ^-constraint 
and (j-constraint respectively, based on the results in Ref. 6. Figures 2 and 3 compare the 
the analytic solutions with the numerical solutions obtained in Ref. 6. We see that both 
the altitude and flight path angle histories match quite well. 
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The trajectory optimization for the aerospace plane stands as a very challenging prob- 
lem because of the highly data-driven model and stringent flight path constraints. Since the 
dominant portion of the optimal ascent trajectory is on the constraints (36) and (37) 3-6 , 
formulas (42)-(45) and their counterparts for the constraint (37) provide a simple yet ac- 
curate representation of 60% — 80% of the trajectory. The trajectory optimization problem 
can be considerably simplified by only numerically investigating the rest of the trajectory. 

It is also worth pointing out that Eq. (32) indicates a nonzero flight path angle along 
the constraints (36) and (37). Applications of the standard energy state approximation 
and zeroth order time-scale decomposition technique to the unconstrained trajectory lead 
to the reduced solution 7 = 0 3,5 . The presence of the two constraints obviously has altered 
the structure of the reduced solution, as suggested by Moerder et al 5 from a different 
perspective. Equations (18) and (19) imply that on the constraint boundary only the 
flight path angle dynamics should be considered “fast” instead of both altitude and flight 
path angle dynamics being considered “fast”. This phenomenon is more evident in the 
next problem. 


3.2 Hypersonic Cruise 

During the early flight tests of the aerospace plane and some other conceivable mis- 
sions, a considerable portion of the trajectory will be hypersonic cruise inside the atmosphere 8 . 
A natural design of the cruising trajectory is a nearly horizontal path with lift-to-drag ra- 
tio nearly maximized. To this end, we would want the cruising trajectory to follow a 
constraint 


Por 0 SC L _ 2 1 / 

C(v ’ p ’ h)= -i^r pv +( 


1 


1 + h (1 + h) 2 


) = 0 


(46) 


where S is the reference area of the vehicle and the lift coefficient corresponding to a*, 
the angle of attack at which the lift-to-drag ratio L/D is maximized. Constraint (46) is 
the result of letting 7 = 0 and L = L* in the right hand side of the 7 ' equation (5) and 
setting 7 ' to zero. It should be noted that during the flight the actual 7 will not be exactly 
zero since the velocity is varying thus o will be slightly different from a*. However, the 
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flight path angle 7 will be near zero and a will be near a* if (46) is followed. Solving v in 
terms of p and h produces 


v = 


(1 + h)( 1 + G(1 + h)p ) 


with the constant G being defined as 

Q _ P^oSCj 

2 m 0 


(47) 


(48) 


Because of the constraint, 7 is expected to be small, and so is h. Therefore, 

(49) 

We first investigate the trajectory numerically. To fly the vehicle, the throttle setting 
is fixed at a constant by specifying a fuel equivalence ratio. The aerodynamic control, a, 
required to follow the constraint (47) is obtained by solving the equation C' = 0. Then 
the a is used in Eqs.(4) and (5) to generate the actual flight trajectory. For computation, 
we have chosen the following data: 


equivalence ratio = 0.5, initial altitude = 40 km, m a = 100,000 kg, 7 0 — 0 


The initial Mach number is 14.63 determined by (47). Figure 4 shows the variations of 7 
and a for the first 10 seconds of the flight obtained through numerical integration. Both 
histories exhibit an unmistakable transient period of 3-4 seconds. After this period 7 and 
a settle down in a steady state with a approaching a* which is 3.93° for our case. This 
asymptotic (as t — > 00) behavior is found to be present regardless of the starting altitude 
and throttle setting, which renders strong support to the theory developed in Section 2.2, 
Eq. (19) in particular. Then we turn to the analytical asymptotic solution of the trajectory 
which is given by Eqs. (35) and (36) with A = l, B = G, n — —1/2 and m = 1: 


7 = 2P( 


e 2z(r) - Q 

e 2z(r) + Q 


) 2 


h = ln( 


G(e 2z(r) +Q ) 2 

—4Qe 2z G) ' 


(50) 

(51) 
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The function z(t) 
period, 


can be calculated as follows. Since a « a* except for the short initial 


V 


K 


D ~ { D )maz 

T C T 
~~ D * SC* D 


(52) 

(53) 


Ct is treated as a constant for the given equivalence ratio and the range of altitude and 
Mach number variations. Then, 


P — = constant, 


z(t ) = Pt 


(54) 


The complete trajectory is designed in the following way: The climbing trajectory 
from takeoff to the beginning of the cruise is numerically found by the inverse dynamics 
technique presented in Ref. 6. This results in a fuel-optimal trajectory subject to inequality 
constraints q < q max and Q < Q max , equality constraints hf = 40 km, 7/ = 0° and 
Vf =Mach 14.63. The cruising trajectory then is analytically modeled by Eqs. (24), (50) 
ari d (51). Figures 5 demonstrates the complete trajectory. The asymptotic solution shows 
an excellent agreement with the cruising trajectory obtained by numerical integration. 
Figure 6 compares the asymptotic and numerical solutions for 7 along the cruising part of 
the trajectory. Except for the a few initial seconds, the asymptotic solution again is very 
close the the true solution. The Mach number increases from 14.63 to 15.29 at the end 
of the ten-minute cruise. It should be noted that the choice of the cruising altitude and 
throttle setting only influences the magnitudes of the h, v and 7. The characteristics of 
the trajectory remain the same. 


Conclusions 

The flight trajectory of aerospace vehicles subject to a class of path constraints has 
been studied. The analysis reveals that under some fairly general conditions the altitude 
dynamics and flight path angle dynamics constitute a natural two-time- scale system: the 
flight path angle dynamics is fast and the altitude dynamics slow. The approximate asymp- 
totic solution for the flight path angle is given as a function of the altitude from which the 
velocity can be expressed as an explicit function of time, regardless of the specific forms of 
the constraints. If the altitude can be solved in terms of the velocity from the constraint, 
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both the altitude and the flight path angle have analytical expressions as functions of time. 
Applications in trajectory analysis for an aerospace plane axe presented. The challenging 
problem of ascent trajectory optimization for the aerospace plane can be significantly sim- 
plified by using these analytical formulas. The technique is also used to design a hypersonic 
cruising trajectory for the aerospace plane. The results strongly support the theory and 
the analytical solutions are in excellent agreement with the numerical results. 
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